Social Control - Final Models

knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)
library(tidybayes)
library(ggplot2)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)

system("shutdown /a")

## [1] 1116

# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE 

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = paste0('_final_with_plan_and_barriers_INTERACTION', as.character(iterations))

df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan', # todo: do we need this??
    'barriers_self_cw', 
    'facilitators_self_cw', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'barriers_self_cb', 
    'facilitators_self_cb', 
    'weartime_self_cb',
    
    # Interactions
    
    'persuasion_self_cw:barriers_self_cw',
    'barriers_self_cw:pressure_self_cw',
    'barriers_self_cw:pushing_self_cw',
    
    'persuasion_self_cw:facilitators_self_cw',
    'pressure_self_cw:facilitators_self_cw',
    'pushing_self_cw:facilitators_self_cw'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  
  'sd(barriers_self_cw)', 
  'sd(facilitators_self_cw)',
  
  'sd(persuasion_self_cw:barriers_self_cw)',
  'sd(barriers_self_cw:pressure_self_cw)',
  'sd(barriers_self_cw:pushing_self_cw)',
  
  'sd(persuasion_self_cw:facilitators_self_cw)',
  'sd(pressure_self_cw:facilitators_self_cw)',
  'sd(pushing_self_cw:facilitators_self_cw)',
  
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')

# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily individual's experienced persuasion",  
    "Daily partner's experienced persuasion", 
    "Daily individual's experienced pressure", 
    "Daily partner's experienced pressure", 
    "Daily individual's experienced pushing", 
    "Daily partner's experienced pushing", 
    "Day", 
    "Own actionplan",
    'Partner actionplan',
    "Daily barriers",
    "Daily facilitators",
    "Daily wear time",
    

    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean individual's experienced persuasion", 
    "Mean partner's experienced persuasion", 
    "Mean individual's experienced pressure", 
    "Mean partner's experienced pressure", 
    "Mean individual's experienced pushing", 
    "Mean partner's experienced pushing", 
    'Mean barriers',
    "Mean facilitators",
    "Mean wear time",
    
    # Interactions
    'persuasion_self_cw:barriers_self_cw',
    'barriers_self_cw:pressure_self_cw',
    'barriers_self_cw:pushing_self_cw',
    
    'persuasion_self_cw:facilitators_self_cw',
    'pressure_self_cw:facilitators_self_cw',
    'pushing_self_cw:facilitators_self_cw'
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily individual's experienced persuasion)", 
  "sd(Daily partner's experienced persuasion)", # OR partner received
  "sd(Daily individual's experienced pressure)", 
  "sd(Daily partner's experienced pressure)", 
  "sd(Daily individual's experienced pushing)", 
  "sd(Daily partner's experienced pushing)", 
  
  'sd(barriers_self_cw)', 
  'sd(facilitators_self_cw)',
  
  'sd(persuasion_self_cw:barriers_self_cw)',
  'sd(barriers_self_cw:pressure_self_cw)',
  'sd(barriers_self_cw:pushing_self_cw)',
  
  'sd(persuasion_self_cw:facilitators_self_cw)',
  'sd(pressure_self_cw:facilitators_self_cw)',
  'sd(pushing_self_cw:facilitators_self_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')

rows_to_pack <- list(
  "Within-Person Effects" = c(2,13),
  "Between-Person Effects" = c(14,22),
  "Interactions" = c(23, 28),
  "Random Effects" = c(29, 43), 
  "Additional Parameters" = c(44,44)
  )


rows_to_pack_ordinal <- list(
  "Within-Person Effects" = c(2+5,13+5),
  "Between-Person Effects" = c(14+5,22+5),
  "Interactions" = c(23+5, 28+5),
  "Random Effects" = c(29+5, 43+5), 
  "Additional Parameters" = c(44+5,44+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 

## [1]   0 720

hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day +  
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | dd | coupleID),
  
  hu = ~ persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | dd | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  control = list(adapt_delta = 0.90),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

pa_sub_digest <- digest::digest(pa_sub)

if (check_models) {
  check_brms(pa_sub, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
}

## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.30 [1.29, 1.30]         1.14      0.77
##       barriers_self_cw 1.15 [1.14, 1.15]         1.07      0.87
##       pressure_self_cw 1.08 [1.08, 1.09]         1.04      0.92
##        pushing_self_cw 1.09 [1.09, 1.09]         1.04      0.92
##   facilitators_self_cw 1.12 [1.12, 1.13]         1.06      0.89
##  persuasion_partner_cw 1.11 [1.10, 1.11]         1.05      0.90
##    pressure_partner_cw 1.05 [1.05, 1.05]         1.02      0.95
##     pushing_partner_cw 1.04 [1.04, 1.05]         1.02      0.96
##     persuasion_self_cb 1.05 [1.05, 1.05]         1.02      0.95
##  persuasion_partner_cb 3.68 [3.66, 3.70]         1.92      0.27
##       pressure_self_cb 3.60 [3.58, 3.62]         1.90      0.28
##    pressure_partner_cb 2.20 [2.19, 2.21]         1.48      0.45
##        pushing_self_cb 2.15 [2.14, 2.16]         1.47      0.47
##     pushing_partner_cb 3.01 [2.99, 3.02]         1.73      0.33
##       barriers_self_cb 3.05 [3.03, 3.07]         1.75      0.33
##   facilitators_self_cb 1.13 [1.12, 1.13]         1.06      0.89
##                    day 1.07 [1.07, 1.08]         1.04      0.93
##  Tolerance 95% CI
##      [0.77, 0.77]
##      [0.87, 0.88]
##      [0.92, 0.93]
##      [0.91, 0.92]
##      [0.89, 0.89]
##      [0.90, 0.91]
##      [0.95, 0.96]
##      [0.95, 0.96]
##      [0.95, 0.96]
##      [0.27, 0.27]
##      [0.28, 0.28]
##      [0.45, 0.46]
##      [0.46, 0.47]
##      [0.33, 0.33]
##      [0.33, 0.33]
##      [0.88, 0.89]
##      [0.93, 0.94]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         94%
##     lognormal          6%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         84%
##               beta-binomial          9%
##                   lognormal          6%
## 
## Divergences:

## 0 of 40000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 40000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 47 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 6, observations = 1776, p-value =
## 0.14
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.003378378
## sample estimates:
## outlier frequency (expected: 0.00164414414414414 ) 
##                                        0.003378378

if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}

# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 

## Sampling priors, please wait...

## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report,
## rope_range = rope_range_continuous, : Coefficients were exponentiated.
## Double check if this was intended.

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack)

	exp(Est.)_hu	SE_hu	95% CI_hu	pd_hu	ROPE_hu	inside ROPE_hu	BF_hu	BF_Evidence_hu	Rhat_hu	Bulk_ESS_hu	Tail_ESS_hu	exp(Est.)_nonzero	SE_nonzero	95% CI_nonzero	pd_nonzero	ROPE_nonzero	inside ROPE_nonzero	BF_nonzero	BF_Evidence_nonzero	Rhat_nonzero	Bulk_ESS_nonzero	Tail_ESS_nonzero
Intercept	3.41***	1.09	[ 1.79, 6.47]	1.000	[0.84, 1.20]	0.001	55.719	Very Strong Evidence	1.000	25963	28327	50.01***	4.20	[42.13, 59.09]	1.000	[0.93, 1.08]	0.000	>100	Overwhelming Evidence	1.000	8663	15562
Within-Person Effects
Daily individual’s experienced persuasion	1.87***	0.28	[ 1.41, 2.61]	1.000	[0.84, 1.20]	0.001	>100	Overwhelming Evidence	1.000	32225	28628	1.00	0.03	[ 0.95, 1.06]	0.569	[0.93, 1.08]	0.993	0.010	Very Strong Evidence for Null	1.000	26933	31420
Daily partner’s experienced persuasion	1.40**	0.17	[ 1.12, 1.83]	0.998	[0.84, 1.20]	0.085	5.292	Moderate Evidence	1.000	32648	26113	1.02	0.02	[ 0.98, 1.07]	0.827	[0.93, 1.08]	0.992	0.013	Very Strong Evidence for Null	1.000	35756	32079
Daily individual’s experienced pressure	0.86	0.33	[ 0.39, 2.09]	0.653	[0.84, 1.20]	0.328	0.213	Moderate Evidence for Null	1.000	31627	24535	0.91	0.06	[ 0.78, 1.05]	0.911	[0.93, 1.08]	0.376	0.028	Very Strong Evidence for Null	1.000	30743	28087
Daily partner’s experienced pressure	1.69	0.72	[ 0.85, 5.72]	0.925	[0.84, 1.20]	0.162	0.451	Weak Evidence for Null	1.000	23866	22260	0.97	0.05	[ 0.88, 1.06]	0.757	[0.93, 1.08]	0.814	0.010	Very Strong Evidence for Null	1.000	40054	34187
Daily individual’s experienced pushing	1.24	0.28	[ 0.77, 2.01]	0.828	[0.84, 1.20]	0.389	0.188	Moderate Evidence for Null	1.000	29899	28023	0.99	0.03	[ 0.93, 1.06]	0.599	[0.93, 1.08]	0.972	0.009	Very Strong Evidence for Null	1.000	35123	35962
Daily partner’s experienced pushing	1.60**	0.30	[ 1.13, 2.54]	0.996	[0.84, 1.20]	0.048	3.367	Moderate Evidence	1.000	28148	22210	0.96	0.03	[ 0.90, 1.02]	0.911	[0.93, 1.08]	0.895	0.020	Very Strong Evidence for Null	1.000	37908	32301
Day	0.91	0.24	[ 0.55, 1.51]	0.636	[0.84, 1.20]	0.480	0.140	Moderate Evidence for Null	1.000	57086	30796	0.95	0.07	[ 0.83, 1.09]	0.745	[0.93, 1.08]	0.630	0.010	Very Strong Evidence for Null	1.000	54096	31338
Own actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Partner actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily barriers	1.14	0.33	[ 0.64, 2.04]	0.670	[0.84, 1.20]	0.434	0.160	Moderate Evidence for Null	1.000	25852	28897	1.26***	0.07	[ 1.13, 1.41]	1.000	[0.93, 1.08]	0.003	20.432	Strong Evidence	1.000	22243	28376
Daily facilitators	1.17	0.31	[ 0.68, 1.96]	0.724	[0.84, 1.20]	0.433	0.158	Moderate Evidence for Null	1.000	26886	28030	1.19**	0.07	[ 1.06, 1.33]	0.997	[0.93, 1.08]	0.054	0.970	Weak Evidence for Null	1.000	17913	23999
Daily wear time	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean individual’s experienced persuasion	1.07	0.73	[ 0.27, 4.02]	0.538	[0.84, 1.20]	0.207	0.338	Weak Evidence for Null	1.000	18162	26070	1.06	0.18	[ 0.75, 1.48]	0.633	[0.93, 1.08]	0.333	0.017	Very Strong Evidence for Null	1.000	16301	24978
Mean partner’s experienced persuasion	1.43	0.95	[ 0.38, 5.24]	0.706	[0.84, 1.20]	0.183	0.388	Weak Evidence for Null	1.000	18874	24215	0.96	0.16	[ 0.69, 1.32]	0.608	[0.93, 1.08]	0.347	0.016	Very Strong Evidence for Null	1.000	15922	24197
Mean individual’s experienced pressure	0.18	0.17	[ 0.03, 1.11]	0.969	[0.84, 1.20]	0.028	2.555	Weak Evidence	1.000	27633	30417	0.68	0.28	[ 0.29, 1.58]	0.823	[0.93, 1.08]	0.091	0.039	Strong Evidence for Null	1.001	6604	11985
Mean partner’s experienced pressure	0.19	0.18	[ 0.03, 1.17]	0.964	[0.84, 1.20]	0.031	2.363	Weak Evidence	1.000	25122	29791	0.49	0.21	[ 0.21, 1.15]	0.950	[0.93, 1.08]	0.035	0.102	Moderate Evidence for Null	1.001	6296	10725
Mean individual’s experienced pushing	1.32	1.26	[ 0.20, 8.70]	0.613	[0.84, 1.20]	0.143	0.485	Weak Evidence for Null	1.000	25432	29055	1.53	0.46	[ 0.84, 2.80]	0.920	[0.93, 1.08]	0.075	0.047	Strong Evidence for Null	1.000	8751	16967
Mean partner’s experienced pushing	1.45	1.41	[ 0.22, 9.51]	0.650	[0.84, 1.20]	0.137	0.517	Weak Evidence for Null	1.000	25870	29954	1.65	0.50	[ 0.91, 3.01]	0.950	[0.93, 1.08]	0.052	0.066	Strong Evidence for Null	1.000	8285	15361
Mean barriers	0.50**	0.13	[ 0.30, 0.83]	0.997	[0.84, 1.20]	0.023	5.207	Moderate Evidence	1.000	56783	30659	0.97	0.07	[ 0.83, 1.12]	0.674	[0.93, 1.08]	0.646	0.011	Very Strong Evidence for Null	1.000	40862	32230
Mean facilitators	1.24	0.23	[ 0.86, 1.80]	0.872	[0.84, 1.20]	0.415	0.175	Moderate Evidence for Null	1.000	54386	31332	1.08	0.05	[ 0.98, 1.19]	0.948	[0.93, 1.08]	0.479	0.049	Strong Evidence for Null	1.000	27594	29862
Mean wear time	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Interactions
persuasion_self_cw:barriers_self_cw	1.02	0.22	[ 0.65, 1.63]	0.528	[0.84, 1.20]	0.581	0.111	Moderate Evidence for Null	1.000	32162	28647	1.03	0.04	[ 0.96, 1.11]	0.777	[0.93, 1.08]	0.903	0.012	Very Strong Evidence for Null	1.000	38494	32777
barriers_self_cw:pressure_self_cw	0.34	0.22	[ 0.09, 1.21]	0.953	[0.84, 1.20]	0.054	1.436	Weak Evidence	1.000	34267	31019	0.95	0.11	[ 0.75, 1.20]	0.684	[0.93, 1.08]	0.441	0.012	Very Strong Evidence for Null	1.000	32615	31479
barriers_self_cw:pushing_self_cw	2.50*	1.01	[ 1.14, 6.03]	0.988	[0.84, 1.20]	0.029	2.808	Weak Evidence	1.000	26846	26558	0.94	0.05	[ 0.85, 1.05]	0.850	[0.93, 1.08]	0.632	0.018	Very Strong Evidence for Null	1.000	30501	30647
persuasion_self_cw:facilitators_self_cw	1.29	0.27	[ 0.87, 2.10]	0.890	[0.84, 1.20]	0.351	0.215	Moderate Evidence for Null	1.000	28751	26602	1.02	0.03	[ 0.96, 1.09]	0.735	[0.93, 1.08]	0.957	0.012	Very Strong Evidence for Null	1.000	27793	27637
pressure_self_cw:facilitators_self_cw	0.88	0.67	[ 0.18, 4.87]	0.568	[0.84, 1.20]	0.185	0.376	Weak Evidence for Null	1.000	29120	26540	0.94	0.10	[ 0.73, 1.18]	0.711	[0.93, 1.08]	0.462	0.014	Very Strong Evidence for Null	1.000	23694	20075
pushing_self_cw:facilitators_self_cw	1.54	0.57	[ 0.77, 3.92]	0.898	[0.84, 1.20]	0.190	0.383	Weak Evidence for Null	1.000	22894	22759	0.95	0.04	[ 0.88, 1.04]	0.879	[0.93, 1.08]	0.757	0.019	Very Strong Evidence for Null	1.000	32770	28974
Random Effects
sd(Intercept)	1.64	0.25	[1.23, 2.23]	NA	NA	NA	NA	NA	1.000	23312	27758	0.32	0.05	[0.25, 0.43]	NA	NA	NA	NA	NA	1.000	14324	22784
sd(Daily individual’s experienced persuasion)	0.21	0.17	[0.01, 0.61]	NA	NA	NA	NA	NA	1.000	15715	20575	0.08	0.03	[0.03, 0.14]	NA	NA	NA	NA	NA	1.000	9028	6570
sd(Daily partner’s experienced persuasion)	0.30	0.17	[0.03, 0.66]	NA	NA	NA	NA	NA	1.000	12444	13969	0.06	0.02	[0.01, 0.11]	NA	NA	NA	NA	NA	1.000	12117	8143
sd(Daily individual’s experienced pressure)	0.54	0.49	[0.02, 2.15]	NA	NA	NA	NA	NA	1.000	17156	20291	0.08	0.08	[0.00, 0.33]	NA	NA	NA	NA	NA	1.001	16273	21030
sd(Daily partner’s experienced pressure)	0.77	0.59	[0.04, 2.37]	NA	NA	NA	NA	NA	1.000	13784	17928	0.05	0.05	[0.00, 0.18]	NA	NA	NA	NA	NA	1.000	20790	19556
sd(Daily individual’s experienced pushing)	0.70	0.36	[0.07, 1.53]	NA	NA	NA	NA	NA	1.000	9718	10941	0.05	0.04	[0.00, 0.14]	NA	NA	NA	NA	NA	1.000	13373	17490
sd(Daily partner’s experienced pushing)	0.44	0.31	[0.03, 1.23]	NA	NA	NA	NA	NA	1.000	11430	16990	0.06	0.04	[0.00, 0.14]	NA	NA	NA	NA	NA	1.000	11705	14178
sd(barriers_self_cw)	1.51	0.24	[1.10, 2.08]	NA	NA	NA	NA	NA	1.000	30311	32165	0.23	0.05	[0.14, 0.35]	NA	NA	NA	NA	NA	1.000	14821	17980
sd(facilitators_self_cw)	1.30	0.27	[0.84, 1.95]	NA	NA	NA	NA	NA	1.000	24583	30170	0.27	0.05	[0.19, 0.39]	NA	NA	NA	NA	NA	1.000	18417	26891
sd(persuasion_self_cw:barriers_self_cw)	0.58	0.33	[0.05, 1.26]	NA	NA	NA	NA	NA	1.000	13106	13548	0.07	0.05	[0.00, 0.18]	NA	NA	NA	NA	NA	1.001	9437	17501
sd(barriers_self_cw:pressure_self_cw)	0.67	0.61	[0.03, 2.58]	NA	NA	NA	NA	NA	1.000	20813	20741	0.11	0.10	[0.01, 0.42]	NA	NA	NA	NA	NA	1.000	21709	23743
sd(barriers_self_cw:pushing_self_cw)	1.26	0.49	[0.38, 2.39]	NA	NA	NA	NA	NA	1.000	14916	14758	0.15	0.05	[0.05, 0.26]	NA	NA	NA	NA	NA	1.000	14581	10756
sd(persuasion_self_cw:facilitators_self_cw)	0.59	0.27	[0.11, 1.23]	NA	NA	NA	NA	NA	1.000	13497	10734	0.07	0.04	[0.00, 0.16]	NA	NA	NA	NA	NA	1.000	9717	13565
sd(pressure_self_cw:facilitators_self_cw)	1.60	1.16	[0.10, 4.28]	NA	NA	NA	NA	NA	1.000	13350	18455	0.12	0.11	[0.01, 0.51]	NA	NA	NA	NA	NA	1.000	16895	22066
sd(pushing_self_cw:facilitators_self_cw)	0.94	0.62	[0.08, 2.45]	NA	NA	NA	NA	NA	1.001	9901	14021	0.05	0.04	[0.00, 0.17]	NA	NA	NA	NA	NA	1.000	19197	21830
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.62	0.01	[0.59, 0.65]	NA	NA	NA	NA	NA	1.000	37840	29108

# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()

## Possible multicollinearity between b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.77), b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.75), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.83), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.71). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

Hurdle part of the model on the left, non-zero part towards the right side of the table

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)

## Scanning ttf files in C:\Windows/Fonts, C:\Users\pascku\AppData\Local/Microsoft/Windows/Fonts ...

## Extracting .afm files from .ttf files...

## C:\Windows\Fonts\Candara.ttf : Candara already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarab.ttf : Candara-Bold already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarai.ttf : Candara-Italic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaral.ttf : Candara-Light already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarali.ttf : Candara-LightItalic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaraz.ttf : Candara-BoldItalic already registered in fonts database. Skipping.
## Found FontName for 0 fonts.
## Scanning afm files in C:/Users/pascku/AppData/Local/R/cache/R/renv/cache/v5/windows/R-4.4/x86_64-w64-mingw32/extrafontdb/1.0/a861555ddec7451c653b40e713166c6f/extrafontdb/metrics

## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.

print(conds_eff)

$persuasion_self_cw

## Warning: Removed 98 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Warning: Removed 191 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Picking joint bandwidth of 0.013

$persuasion_partner_cw

## Warning: Removed 118 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Warning: Removed 240 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Picking joint bandwidth of 0.00825

$pressure_self_cw

## Warning: Removed 76 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Warning: Removed 117 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Picking joint bandwidth of 0.0185

$pressure_partner_cw

## Warning: Removed 75 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Warning: Removed 100 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Picking joint bandwidth of 0.0326

$pushing_self_cw

## Warning: Removed 100 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Warning: Removed 58 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Picking joint bandwidth of 0.0136

$pushing_partner_cw

## Warning: Removed 35 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Warning: Removed 75 rows containing missing values or values outside the scale
## range (`geom_line()`).

## Picking joint bandwidth of 0.014

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

Comparing effect size of pressure and pushing

hypothesis(pa_sub, "pressure_self_cw < pushing_self_cw")

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... < 0    -0.09      0.09    -0.23     0.05       5.71
##   Post.Prob Star
## 1      0.85     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 

## [1]   5.75 971.25

hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  control = list(adapt_delta = 0.99),
  iter = iterations + 200,
  warmup = warmup + 200,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

pa_obj_log_digest <- digest::digest(pa_obj_log)

if (check_models) {
  check_brms(pa_obj_log, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
}

## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.20 [1.16, 1.24]         1.09      0.84
##       barriers_self_cw 1.20 [1.17, 1.24]         1.10      0.83
##       pressure_self_cw 1.17 [1.13, 1.21]         1.08      0.86
##        pushing_self_cw 1.26 [1.22, 1.30]         1.12      0.80
##   facilitators_self_cw 1.17 [1.14, 1.21]         1.08      0.85
##  persuasion_partner_cw 1.14 [1.11, 1.18]         1.07      0.88
##    pressure_partner_cw 1.10 [1.08, 1.14]         1.05      0.91
##     pushing_partner_cw 1.16 [1.12, 1.20]         1.08      0.86
##       pressure_self_cb 3.80 [3.62, 3.98]         1.95      0.26
##    pressure_partner_cb 4.01 [3.83, 4.21]         2.00      0.25
##       barriers_self_cb 1.18 [1.15, 1.22]         1.09      0.85
##   facilitators_self_cb 1.16 [1.13, 1.20]         1.08      0.86
##                    day 1.06 [1.04, 1.10]         1.03      0.94
##       weartime_self_cw 1.06 [1.04, 1.10]         1.03      0.94
##       weartime_self_cb 1.21 [1.18, 1.26]         1.10      0.82
##  Tolerance 95% CI
##      [0.81, 0.86]
##      [0.80, 0.86]
##      [0.83, 0.88]
##      [0.77, 0.82]
##      [0.83, 0.88]
##      [0.85, 0.90]
##      [0.88, 0.93]
##      [0.84, 0.89]
##      [0.25, 0.28]
##      [0.24, 0.26]
##      [0.82, 0.87]
##      [0.83, 0.89]
##      [0.91, 0.96]
##      [0.91, 0.96]
##      [0.80, 0.85]
## 
## Moderate Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 7.76 [7.38, 8.17]         2.79      0.13
##  persuasion_partner_cb 8.17 [7.77, 8.60]         2.86      0.12
##        pushing_self_cb 5.88 [5.60, 6.18]         2.42      0.17
##     pushing_partner_cb 5.74 [5.47, 6.03]         2.40      0.17
##  Tolerance 95% CI
##      [0.12, 0.14]
##      [0.12, 0.13]
##      [0.16, 0.18]
##      [0.17, 0.18]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         88%
##     lognormal          9%
##       weibull          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##                   lognormal         72%
##                     tweedie         16%
##  neg. binomial (zero-infl.)          9%
## 
## Divergences:

## 0 of 40000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 40000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 11 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 14, observations = 1594, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0002979925 0.0056461731
## sample estimates:
## outlier frequency (expected: 0.003099121706399 ) 
##                                      0.008782936

if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}

# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 

## Sampling priors, please wait...

## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report, :
## Coefficients were exponentiated. Double check if this was intended.

summary_pa_obj %>%
  print_df(rows_to_pack = rows_to_pack)

	exp(Est.)	SE	95% CI	pd	ROPE	inside ROPE	BF	BF_Evidence	Rhat	Bulk_ESS	Tail_ESS
Intercept	123.06***	7.77	[108.16, 139.27]	1.000	[0.94, 1.06]	0.000	>100	Overwhelming Evidence	1.000	8031	14512
Within-Person Effects
Daily individual’s experienced persuasion	1.02	0.02	[ 0.98, 1.06]	0.864	[0.94, 1.06]	0.985	0.013	Very Strong Evidence for Null	1.000	34459	31842
Daily partner’s experienced persuasion	1.03	0.02	[ 0.99, 1.07]	0.931	[0.94, 1.06]	0.969	0.021	Very Strong Evidence for Null	1.000	38113	30960
Daily individual’s experienced pressure	0.97	0.05	[ 0.87, 1.07]	0.759	[0.94, 1.06]	0.672	0.009	Very Strong Evidence for Null	1.000	40938	32241
Daily partner’s experienced pressure	0.97	0.04	[ 0.90, 1.05]	0.793	[0.94, 1.06]	0.797	0.009	Very Strong Evidence for Null	1.000	48631	31224
Daily individual’s experienced pushing	1.02	0.03	[ 0.96, 1.08]	0.734	[0.94, 1.06]	0.939	0.009	Very Strong Evidence for Null	1.000	37349	33850
Daily partner’s experienced pushing	1.00	0.02	[ 0.96, 1.05]	0.557	[0.94, 1.06]	0.988	0.006	Very Strong Evidence for Null	1.000	43220	32281
Day	0.96	0.05	[ 0.87, 1.06]	0.771	[0.94, 1.06]	0.666	0.008	Very Strong Evidence for Null	1.000	60635	30163
Own actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Partner actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily barriers	1.01	0.03	[ 0.94, 1.07]	0.562	[0.94, 1.06]	0.939	0.008	Very Strong Evidence for Null	1.000	28586	29795
Daily facilitators	1.04	0.03	[ 0.99, 1.10]	0.945	[0.94, 1.06]	0.778	0.031	Strong Evidence for Null	1.000	31594	27535
Daily wear time	1.00***	0.00	[ 1.00, 1.00]	1.000	[0.94, 1.06]	1.000	2.621	Weak Evidence	1.000	66527	31035
Between-Person Effects
Mean individual’s experienced persuasion	0.99	0.17	[ 0.69, 1.40]	0.533	[0.94, 1.06]	0.274	0.016	Very Strong Evidence for Null	1.001	6637	13062
Mean partner’s experienced persuasion	0.86	0.15	[ 0.61, 1.23]	0.796	[0.94, 1.06]	0.196	0.023	Very Strong Evidence for Null	1.001	6727	12797
Mean individual’s experienced pressure	1.12	0.21	[ 0.77, 1.64]	0.730	[0.94, 1.06]	0.216	0.012	Very Strong Evidence for Null	1.001	9668	17120
Mean partner’s experienced pressure	1.04	0.19	[ 0.72, 1.49]	0.579	[0.94, 1.06]	0.261	0.011	Very Strong Evidence for Null	1.001	8766	16544
Mean individual’s experienced pushing	1.11	0.27	[ 0.69, 1.82]	0.668	[0.94, 1.06]	0.186	0.014	Very Strong Evidence for Null	1.001	9898	17862
Mean partner’s experienced pushing	1.27	0.31	[ 0.78, 2.05]	0.836	[0.94, 1.06]	0.127	0.022	Very Strong Evidence for Null	1.001	9438	16390
Mean barriers	0.80***	0.04	[ 0.72, 0.89]	1.000	[0.94, 1.06]	0.001	26.021	Strong Evidence	1.000	52877	32143
Mean facilitators	0.96	0.04	[ 0.89, 1.03]	0.873	[0.94, 1.06]	0.704	0.019	Very Strong Evidence for Null	1.000	34078	32619
Mean wear time	1.00*	0.00	[ 1.00, 1.00]	0.988	[0.94, 1.06]	1.000	0.124	Moderate Evidence for Null	1.000	31806	31751
Interactions
persuasion_self_cw:barriers_self_cw	0.97	0.02	[ 0.93, 1.02]	0.882	[0.94, 1.06]	0.905	0.013	Very Strong Evidence for Null	1.000	53504	32603
barriers_self_cw:pressure_self_cw	0.96	0.08	[ 0.81, 1.14]	0.684	[0.94, 1.06]	0.495	0.008	Very Strong Evidence for Null	1.000	34654	33293
barriers_self_cw:pushing_self_cw	1.02	0.03	[ 0.95, 1.09]	0.696	[0.94, 1.06]	0.894	0.007	Very Strong Evidence for Null	1.000	43986	31833
persuasion_self_cw:facilitators_self_cw	0.99	0.02	[ 0.95, 1.04]	0.629	[0.94, 1.06]	0.994	0.008	Very Strong Evidence for Null	1.000	38120	29806
pressure_self_cw:facilitators_self_cw	0.92	0.10	[ 0.69, 1.14]	0.798	[0.94, 1.06]	0.328	0.016	Very Strong Evidence for Null	1.000	19016	20911
pushing_self_cw:facilitators_self_cw	1.03	0.03	[ 0.96, 1.09]	0.793	[0.94, 1.06]	0.868	0.011	Very Strong Evidence for Null	1.000	36408	28457
Random Effects
sd(Intercept)	0.34	0.05	[0.26, 0.45]	NA	NA	NA	NA	NA	1.000	10470	18026
sd(Daily individual’s experienced persuasion)	0.05	0.02	[0.01, 0.09]	NA	NA	NA	NA	NA	1.000	12921	10061
sd(Daily partner’s experienced persuasion)	0.05	0.02	[0.01, 0.10]	NA	NA	NA	NA	NA	1.000	10508	9922
sd(Daily individual’s experienced pressure)	0.05	0.05	[0.00, 0.20]	NA	NA	NA	NA	NA	1.000	16988	19408
sd(Daily partner’s experienced pressure)	0.03	0.03	[0.00, 0.12]	NA	NA	NA	NA	NA	1.000	22902	19193
sd(Daily individual’s experienced pushing)	0.08	0.04	[0.01, 0.15]	NA	NA	NA	NA	NA	1.000	10344	11953
sd(Daily partner’s experienced pushing)	0.04	0.03	[0.00, 0.11]	NA	NA	NA	NA	NA	1.000	12602	15278
sd(barriers_self_cw)	0.13	0.03	[0.06, 0.20]	NA	NA	NA	NA	NA	1.000	11395	8797
sd(facilitators_self_cw)	0.07	0.04	[0.01, 0.15]	NA	NA	NA	NA	NA	1.000	9088	10195
sd(persuasion_self_cw:barriers_self_cw)	0.03	0.02	[0.00, 0.09]	NA	NA	NA	NA	NA	1.000	19050	20948
sd(barriers_self_cw:pressure_self_cw)	0.07	0.07	[0.00, 0.31]	NA	NA	NA	NA	NA	1.000	20796	21799
sd(barriers_self_cw:pushing_self_cw)	0.03	0.03	[0.00, 0.10]	NA	NA	NA	NA	NA	1.000	21260	21607
sd(persuasion_self_cw:facilitators_self_cw)	0.03	0.02	[0.00, 0.09]	NA	NA	NA	NA	NA	1.000	13572	18030
sd(pressure_self_cw:facilitators_self_cw)	0.24	0.16	[0.02, 0.62]	NA	NA	NA	NA	NA	1.000	10644	14130
sd(pushing_self_cw:facilitators_self_cw)	0.04	0.04	[0.00, 0.14]	NA	NA	NA	NA	NA	1.000	15627	18450
Additional Parameters
sigma	0.53	0.01	[0.51, 0.55]	NA	NA	NA	NA	NA	1.000	40490	31103

plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()

## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.87), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.72), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.78). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)

## [1] 0 5

hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

mood_gauss_digest <- digest::digest(mood_gauss)

if (check_models) {
  check_brms(mood_gauss, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss, integer = FALSE)
}

## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.24 [1.20,  1.29]         1.11      0.81
##       barriers_self_cw 1.12 [1.09,  1.16]         1.06      0.89
##       pressure_self_cw 1.14 [1.11,  1.18]         1.07      0.87
##        pushing_self_cw 1.30 [1.26,  1.34]         1.14      0.77
##   facilitators_self_cw 1.09 [1.06,  1.13]         1.04      0.92
##  persuasion_partner_cw 1.14 [1.11,  1.19]         1.07      0.87
##    pressure_partner_cw 1.09 [1.06,  1.13]         1.04      0.92
##     pushing_partner_cw 1.13 [1.10,  1.17]         1.06      0.88
##       pressure_self_cb 5.00 [4.75,  5.25]         2.24      0.20
##    pressure_partner_cb 4.87 [4.63,  5.12]         2.21      0.21
##       barriers_self_cb 1.12 [1.09,  1.16]         1.06      0.89
##   facilitators_self_cb 1.09 [1.06,  1.13]         1.04      0.92
##                    day 1.08 [1.06,  1.12]         1.04      0.92
##  Tolerance 95% CI
##      [0.78, 0.83]
##      [0.86, 0.92]
##      [0.84, 0.90]
##      [0.74, 0.80]
##      [0.89, 0.94]
##      [0.84, 0.90]
##      [0.89, 0.94]
##      [0.85, 0.91]
##      [0.19, 0.21]
##      [0.20, 0.22]
##      [0.86, 0.92]
##      [0.88, 0.94]
##      [0.89, 0.95]
## 
## Moderate Correlation
## 
##                Term  VIF    VIF 95% CI Increased SE Tolerance
##     pushing_self_cb 7.80 [7.41,  8.22]         2.79      0.13
##  pushing_partner_cb 7.80 [7.41,  8.22]         2.79      0.13
##  Tolerance 95% CI
##      [0.12, 0.13]
##      [0.12, 0.13]
## 
## High Correlation
## 
##                   Term   VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 10.41 [9.88, 10.98]         3.23      0.10
##  persuasion_partner_cb 10.26 [9.73, 10.81]         3.20      0.10
##  Tolerance 95% CI
##      [0.09, 0.10]
##      [0.09, 0.10]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         41%
##        normal         41%
##   exponential          6%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##        tweedie         41%
##  beta-binomial         38%
##    half-cauchy          6%
## 
## Divergences:

## 0 of 40000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 40000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 10 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 9, observations = 1776, p-value =
## 0.01064
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.002319754 0.009597934
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.005067568

if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}

summary_mood <- summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) 

## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

## Warning in compute_rope(range = rope_range): Collinearity detected. Some
## VIFs are > 10. This may invalidate ROPE inferences!

Check for Multicollinearity

Low Correlation

              Term  VIF    VIF 95% CI Increased SE Tolerance
persuasion_self_cw 1.24 [1.20,  1.29]         1.11      0.81
  barriers_self_cw 1.12 [1.09,  1.16]         1.06      0.89
  pressure_self_cw 1.14 [1.11,  1.18]         1.07      0.87
   pushing_self_cw 1.30 [1.26,  1.34]         1.14      0.77

facilitators_self_cw 1.09 [1.06, 1.13] 1.04 0.92 persuasion_partner_cw 1.14 [1.11, 1.19] 1.07 0.87 pressure_partner_cw 1.09 [1.06, 1.13] 1.04 0.92 pushing_partner_cw 1.13 [1.10, 1.17] 1.06 0.88 pressure_self_cb 5.00 [4.75, 5.25] 2.24 0.20 pressure_partner_cb 4.87 [4.63, 5.12] 2.21 0.21 barriers_self_cb 1.12 [1.09, 1.16] 1.06 0.89 facilitators_self_cb 1.09 [1.06, 1.13] 1.04 0.92 day 1.08 [1.06, 1.12] 1.04 0.92 Tolerance 95% CI [0.78, 0.83] [0.86, 0.92] [0.84, 0.90] [0.74, 0.80] [0.89, 0.94] [0.84, 0.90] [0.89, 0.94] [0.85, 0.91] [0.19, 0.21] [0.20, 0.22] [0.86, 0.92] [0.88, 0.94] [0.89, 0.95]

Moderate Correlation

           Term  VIF    VIF 95% CI Increased SE Tolerance
pushing_self_cb 7.80 [7.41,  8.22]         2.79      0.13

pushing_partner_cb 7.80 [7.41, 8.22] 2.79 0.13 Tolerance 95% CI [0.12, 0.13] [0.12, 0.13]

High Correlation

              Term   VIF    VIF 95% CI Increased SE Tolerance
persuasion_self_cb 10.41 [9.88, 10.98]         3.23      0.10

persuasion_partner_cb 10.26 [9.73, 10.81] 3.20 0.10 Tolerance 95% CI [0.09, 0.10] [0.09, 0.10]

## Sampling priors, please wait...

summary_mood %>%
  print_df(rows_to_pack = rows_to_pack)

	Est.	SE	95% CI	pd	ROPE	inside ROPE	BF	BF_Evidence	Rhat	Bulk_ESS	Tail_ESS
Intercept	3.79***	0.11	[ 3.58, 4.01]	1.000	[-0.11, 0.11]	0.000	>100	Overwhelming Evidence	1.000	8034	15527
Within-Person Effects
Daily individual’s experienced persuasion	0.01	0.03	[-0.04, 0.06]	0.615	[-0.11, 0.11]	1.000	0.005	Very Strong Evidence for Null	1.000	58341	33448
Daily partner’s experienced persuasion	0.01	0.02	[-0.04, 0.06]	0.663	[-0.11, 0.11]	1.000	0.005	Very Strong Evidence for Null	1.000	58314	31804
Daily individual’s experienced pressure	-0.02	0.07	[-0.16, 0.12]	0.619	[-0.11, 0.11]	0.859	0.005	Very Strong Evidence for Null	1.000	47243	34060
Daily partner’s experienced pressure	-0.05	0.06	[-0.18, 0.07]	0.810	[-0.11, 0.11]	0.825	0.008	Very Strong Evidence for Null	1.000	44211	31698
Daily individual’s experienced pushing	0.02	0.04	[-0.06, 0.10]	0.705	[-0.11, 0.11]	0.980	0.006	Very Strong Evidence for Null	1.000	43183	35615
Daily partner’s experienced pushing	0.09*	0.04	[ 0.02, 0.16]	0.990	[-0.11, 0.11]	0.690	0.097	Strong Evidence for Null	1.000	47310	30460
Day	0.26***	0.08	[ 0.11, 0.41]	1.000	[-0.11, 0.11]	0.025	1.134	Weak Evidence	1.000	67179	31237
Own actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Partner actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily barriers	-0.12	0.06	[-0.25, 0.01]	0.969	[-0.11, 0.11]	0.431	0.049	Strong Evidence for Null	1.000	20011	25278
Daily facilitators	0.11**	0.03	[ 0.05, 0.17]	0.999	[-0.11, 0.11]	0.502	1.022	Weak Evidence	1.000	55329	29865
Daily wear time	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean individual’s experienced persuasion	0.13	0.30	[-0.47, 0.72]	0.667	[-0.11, 0.11]	0.265	0.015	Very Strong Evidence for Null	1.001	7278	14121
Mean partner’s experienced persuasion	0.42	0.29	[-0.18, 1.00]	0.918	[-0.11, 0.11]	0.112	0.040	Strong Evidence for Null	1.001	7123	13430
Mean individual’s experienced pressure	-0.07	0.30	[-0.67, 0.53]	0.594	[-0.11, 0.11]	0.275	0.009	Very Strong Evidence for Null	1.000	8986	17404
Mean partner’s experienced pressure	-0.32	0.30	[-0.91, 0.28]	0.854	[-0.11, 0.11]	0.166	0.016	Very Strong Evidence for Null	1.000	8783	16707
Mean individual’s experienced pushing	0.23	0.42	[-0.61, 1.08]	0.705	[-0.11, 0.11]	0.178	0.014	Very Strong Evidence for Null	1.001	10000	16141
Mean partner’s experienced pushing	-0.06	0.42	[-0.90, 0.79]	0.557	[-0.11, 0.11]	0.208	0.012	Very Strong Evidence for Null	1.000	9929	16149
Mean barriers	-0.50***	0.08	[-0.65, -0.34]	1.000	[-0.11, 0.11]	0.000	>100	Overwhelming Evidence	1.000	72710	31917
Mean facilitators	0.29***	0.05	[ 0.18, 0.39]	1.000	[-0.11, 0.11]	0.001	>100	Overwhelming Evidence	1.000	43379	32312
Mean wear time	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Interactions
persuasion_self_cw:barriers_self_cw	0.00	0.04	[-0.08, 0.07]	0.529	[-0.11, 0.11]	0.993	0.005	Very Strong Evidence for Null	1.000	50722	33188
barriers_self_cw:pressure_self_cw	-0.14	0.16	[-0.47, 0.20]	0.808	[-0.11, 0.11]	0.364	0.011	Very Strong Evidence for Null	1.000	31974	31759
barriers_self_cw:pushing_self_cw	0.08	0.05	[-0.02, 0.18]	0.950	[-0.11, 0.11]	0.701	0.019	Very Strong Evidence for Null	1.000	55985	33379
persuasion_self_cw:facilitators_self_cw	0.00	0.03	[-0.06, 0.06]	0.530	[-0.11, 0.11]	0.999	0.005	Very Strong Evidence for Null	1.000	51081	31523
pressure_self_cw:facilitators_self_cw	0.00	0.13	[-0.24, 0.31]	0.502	[-0.11, 0.11]	0.611	0.007	Very Strong Evidence for Null	1.000	27612	25163
pushing_self_cw:facilitators_self_cw	0.02	0.05	[-0.08, 0.13]	0.673	[-0.11, 0.11]	0.935	0.007	Very Strong Evidence for Null	1.000	34592	31350
Random Effects
sd(Intercept)	0.62	0.08	[0.49, 0.81]	NA	NA	NA	NA	NA	1.000	9965	19755
sd(Daily individual’s experienced persuasion)	0.04	0.03	[0.00, 0.11]	NA	NA	NA	NA	NA	1.000	13775	19005
sd(Daily partner’s experienced persuasion)	0.04	0.03	[0.00, 0.11]	NA	NA	NA	NA	NA	1.000	16581	19683
sd(Daily individual’s experienced pressure)	0.07	0.06	[0.00, 0.25]	NA	NA	NA	NA	NA	1.000	23395	22857
sd(Daily partner’s experienced pressure)	0.09	0.07	[0.00, 0.29]	NA	NA	NA	NA	NA	1.000	16353	19990
sd(Daily individual’s experienced pushing)	0.09	0.05	[0.01, 0.20]	NA	NA	NA	NA	NA	1.001	11251	12691
sd(Daily partner’s experienced pushing)	0.09	0.05	[0.01, 0.19]	NA	NA	NA	NA	NA	1.000	13536	16131
sd(barriers_self_cw)	0.32	0.05	[0.23, 0.43]	NA	NA	NA	NA	NA	1.000	19603	27622
sd(facilitators_self_cw)	0.05	0.04	[0.00, 0.16]	NA	NA	NA	NA	NA	1.000	13807	20577
sd(persuasion_self_cw:barriers_self_cw)	0.07	0.05	[0.00, 0.18]	NA	NA	NA	NA	NA	1.000	14309	19001
sd(barriers_self_cw:pressure_self_cw)	0.41	0.18	[0.08, 0.85]	NA	NA	NA	NA	NA	1.000	14368	12763
sd(barriers_self_cw:pushing_self_cw)	0.05	0.05	[0.00, 0.18]	NA	NA	NA	NA	NA	1.000	20198	24464
sd(persuasion_self_cw:facilitators_self_cw)	0.04	0.03	[0.00, 0.13]	NA	NA	NA	NA	NA	1.000	17608	20566
sd(pressure_self_cw:facilitators_self_cw)	0.16	0.13	[0.01, 0.55]	NA	NA	NA	NA	NA	1.000	19179	20271
sd(pushing_self_cw:facilitators_self_cw)	0.11	0.06	[0.01, 0.26]	NA	NA	NA	NA	NA	1.000	14431	15173
Additional Parameters
sigma	0.84	0.02	[0.82, 0.88]	NA	NA	NA	NA	NA	1.000	53758	30528

plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()

## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.76), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.73), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.73), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.75), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.85). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 

## [1] 0 5

hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

reactance_ordinal_digest <- digest::digest(reactance_ordinal)

if (check_models) {
  check_brms(reactance_ordinal)
  DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')
}

## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##       pressure_self_cw 4.08 [3.79,  4.40]         2.02      0.24
##        pushing_self_cw 1.37 [1.31,  1.45]         1.17      0.73
##   facilitators_self_cw 1.27 [1.21,  1.34]         1.13      0.79
##  persuasion_partner_cw 1.69 [1.60,  1.80]         1.30      0.59
##    pressure_partner_cw 1.12 [1.08,  1.19]         1.06      0.89
##     pushing_partner_cw 1.29 [1.23,  1.37]         1.14      0.77
##     persuasion_self_cb 1.91 [1.80,  2.04]         1.38      0.52
##  persuasion_partner_cb 1.22 [1.17,  1.29]         1.11      0.82
##       pressure_self_cb 1.06 [1.03,  1.13]         1.03      0.94
##    pressure_partner_cb 1.23 [1.18,  1.30]         1.11      0.81
##       barriers_self_cb 4.96 [4.60,  5.36]         2.23      0.20
##   facilitators_self_cb 4.96 [4.60,  5.36]         2.23      0.20
##                    day 2.14 [2.01,  2.28]         1.46      0.47
##  Tolerance 95% CI
##      [0.23, 0.26]
##      [0.69, 0.77]
##      [0.75, 0.83]
##      [0.55, 0.62]
##      [0.84, 0.93]
##      [0.73, 0.81]
##      [0.49, 0.56]
##      [0.77, 0.86]
##      [0.88, 0.97]
##      [0.77, 0.85]
##      [0.19, 0.22]
##      [0.19, 0.22]
##      [0.44, 0.50]
## 
## Moderate Correlation
## 
##                Term  VIF    VIF 95% CI Increased SE Tolerance
##  persuasion_self_cw 8.57 [7.91,  9.28]         2.93      0.12
##     pushing_self_cb 5.03 [4.66,  5.43]         2.24      0.20
##  pushing_partner_cb 6.21 [5.75,  6.72]         2.49      0.16
##  Tolerance 95% CI
##      [0.11, 0.13]
##      [0.18, 0.21]
##      [0.15, 0.17]
## 
## High Correlation
## 
##              Term   VIF    VIF 95% CI Increased SE Tolerance
##  barriers_self_cw 10.36 [9.56, 11.24]         3.22      0.10
##  Tolerance 95% CI
##      [0.09, 0.10]
## Error in h(simpleError(msg, call)) : 
##   error in evaluating the argument 'x' in selecting a method for function 'print': Predictive errors are not defined for ordinal or categorical models.
## 
## Divergences:

## 0 of 40000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 40000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 36 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 486, p-value = 0.38
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002057613
## sample estimates:
## outlier frequency (expected: 0.000432098765432099 ) 
##                                         0.002057613

if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}

summary_reactance_ordinal <- summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) 

## Model has interaction terms. VIFs might be inflated.
##   You may check multicollinearity among predictors of a model without
##   interaction terms.

## Warning in compute_rope(range = rope_range): Collinearity detected. Some
## VIFs are > 10. This may invalidate ROPE inferences!

Check for Multicollinearity

Low Correlation

              Term  VIF    VIF 95% CI Increased SE Tolerance
  pressure_self_cw 4.08 [3.79,  4.40]         2.02      0.24
   pushing_self_cw 1.37 [1.31,  1.45]         1.17      0.73

facilitators_self_cw 1.27 [1.21, 1.34] 1.13 0.79 persuasion_partner_cw 1.69 [1.60, 1.80] 1.30 0.59 pressure_partner_cw 1.12 [1.08, 1.19] 1.06 0.89 pushing_partner_cw 1.29 [1.23, 1.37] 1.14 0.77 persuasion_self_cb 1.91 [1.80, 2.04] 1.38 0.52 persuasion_partner_cb 1.22 [1.17, 1.29] 1.11 0.82 pressure_self_cb 1.06 [1.03, 1.13] 1.03 0.94 pressure_partner_cb 1.23 [1.18, 1.30] 1.11 0.81 barriers_self_cb 4.96 [4.60, 5.36] 2.23 0.20 facilitators_self_cb 4.96 [4.60, 5.36] 2.23 0.20 day 2.14 [2.01, 2.28] 1.46 0.47 Tolerance 95% CI [0.23, 0.26] [0.69, 0.77] [0.75, 0.83] [0.55, 0.62] [0.84, 0.93] [0.73, 0.81] [0.49, 0.56] [0.77, 0.86] [0.88, 0.97] [0.77, 0.85] [0.19, 0.22] [0.19, 0.22] [0.44, 0.50]

Moderate Correlation

           Term  VIF    VIF 95% CI Increased SE Tolerance

persuasion_self_cw 8.57 [7.91, 9.28] 2.93 0.12 pushing_self_cb 5.03 [4.66, 5.43] 2.24 0.20 pushing_partner_cb 6.21 [5.75, 6.72] 2.49 0.16 Tolerance 95% CI [0.11, 0.13] [0.18, 0.21] [0.15, 0.17]

High Correlation

         Term   VIF    VIF 95% CI Increased SE Tolerance

barriers_self_cw 10.36 [9.56, 11.24] 3.22 0.10 Tolerance 95% CI [0.09, 0.10]

## Sampling priors, please wait...

summary_reactance_ordinal %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)

	OR	SE	95% CI	pd	ROPE	inside ROPE	BF	BF_Evidence	Rhat	Bulk_ESS	Tail_ESS
Intercept	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Intercept[1]	5.24***	2.33	[ 2.23, 13.19]	1.000	[0.84, 1.20]	0.000	>100	Overwhelming Evidence	1.000	33731	29944
Intercept[2]	13.11***	6.13	[ 5.42, 34.96]	1.000	[0.84, 1.20]	0.000	>100	Overwhelming Evidence	1.000	31987	29285
Intercept[3]	43.61***	22.10	[ 17.01, 128.11]	1.000	[0.84, 1.20]	0.000	>100	Overwhelming Evidence	1.000	29869	29288
Intercept[4]	254.83***	155.72	[ 83.32, 917.86]	1.000	[0.84, 1.20]	0.000	>100	Overwhelming Evidence	1.000	27685	28142
Intercept[5]	13351.51***	16831.35	[1628.97, 285902.89]	1.000	[0.84, 1.20]	0.000	>100	Overwhelming Evidence	1.000	39002	29466
Within-Person Effects
Daily individual’s experienced persuasion	0.64**	0.10	[ 0.46, 0.86]	0.999	[0.84, 1.20]	0.039	7.369	Moderate Evidence	1.000	32737	28343
Daily partner’s experienced persuasion	1.10	0.17	[ 0.80, 1.48]	0.741	[0.84, 1.20]	0.659	0.088	Strong Evidence for Null	1.000	36870	29745
Daily individual’s experienced pressure	2.21*	0.71	[ 1.08, 4.26]	0.983	[0.84, 1.20]	0.037	1.223	Weak Evidence	1.000	26836	26589
Daily partner’s experienced pressure	1.13	0.43	[ 0.41, 2.46]	0.628	[0.84, 1.20]	0.337	0.091	Strong Evidence for Null	1.000	26133	20054
Daily individual’s experienced pushing	1.30	0.20	[ 0.93, 1.78]	0.944	[0.84, 1.20]	0.302	0.271	Moderate Evidence for Null	1.000	38983	28672
Daily partner’s experienced pushing	0.92	0.16	[ 0.65, 1.32]	0.684	[0.84, 1.20]	0.642	0.075	Strong Evidence for Null	1.000	41776	30483
Day	1.33	0.76	[ 0.43, 4.11]	0.688	[0.84, 1.20]	0.217	0.071	Strong Evidence for Null	1.000	45608	32116
Own actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Partner actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily barriers	1.28	0.51	[ 0.58, 2.94]	0.737	[0.84, 1.20]	0.290	0.086	Strong Evidence for Null	1.000	43277	33473
Daily facilitators	1.15	0.51	[ 0.44, 2.67]	0.618	[0.84, 1.20]	0.294	0.093	Strong Evidence for Null	1.000	29957	29903
Daily wear time	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean individual’s experienced persuasion	1.49	1.31	[ 0.25, 9.25]	0.674	[0.84, 1.20]	0.149	0.100	Strong Evidence for Null	1.000	22625	27377
Mean partner’s experienced persuasion	1.94	1.88	[ 0.31, 15.31]	0.762	[0.84, 1.20]	0.120	0.114	Moderate Evidence for Null	1.000	22074	25540
Mean individual’s experienced pressure	2.88	3.01	[ 0.37, 24.66]	0.848	[0.84, 1.20]	0.083	0.128	Moderate Evidence for Null	1.000	26998	30390
Mean partner’s experienced pressure	1.08	1.09	[ 0.13, 6.97]	0.531	[0.84, 1.20]	0.142	0.078	Strong Evidence for Null	1.000	21451	25485
Mean individual’s experienced pushing	1.02	1.29	[ 0.08, 14.46]	0.507	[0.84, 1.20]	0.115	0.086	Strong Evidence for Null	1.000	23211	27282
Mean partner’s experienced pushing	0.02**	0.03	[ 0.00, 0.33]	0.996	[0.84, 1.20]	0.003	3.686	Moderate Evidence	1.000	24552	26788
Mean barriers	0.89	0.63	[ 0.21, 3.48]	0.564	[0.84, 1.20]	0.200	0.081	Strong Evidence for Null	1.000	32718	28853
Mean facilitators	0.60	0.28	[ 0.23, 1.49]	0.871	[0.84, 1.20]	0.168	0.194	Moderate Evidence for Null	1.000	28304	29905
Mean wear time	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Interactions
persuasion_self_cw:barriers_self_cw	0.69	0.15	[ 0.43, 1.04]	0.961	[0.84, 1.20]	0.182	0.365	Weak Evidence for Null	1.000	32999	29006
barriers_self_cw:pressure_self_cw	1.71	0.75	[ 0.70, 4.17]	0.888	[0.84, 1.20]	0.153	0.157	Moderate Evidence for Null	1.000	35588	33208
barriers_self_cw:pushing_self_cw	0.69	0.16	[ 0.42, 1.08]	0.947	[0.84, 1.20]	0.195	0.264	Moderate Evidence for Null	1.000	35035	29498
persuasion_self_cw:facilitators_self_cw	0.91	0.18	[ 0.62, 1.35]	0.682	[0.84, 1.20]	0.586	0.097	Strong Evidence for Null	1.000	33388	30290
pressure_self_cw:facilitators_self_cw	1.68	1.00	[ 0.55, 6.63]	0.820	[0.84, 1.20]	0.168	0.173	Moderate Evidence for Null	1.000	18589	20585
pushing_self_cw:facilitators_self_cw	0.65	0.16	[ 0.40, 1.12]	0.945	[0.84, 1.20]	0.150	0.431	Weak Evidence for Null	1.000	28963	23880
Random Effects
sd(Intercept)	1.30	0.34	[0.74, 2.12]	NA	NA	NA	NA	NA	1.000	13995	22011
sd(Daily individual’s experienced persuasion)	0.34	0.23	[0.02, 0.81]	NA	NA	NA	NA	NA	1.000	8114	15291
sd(Daily partner’s experienced persuasion)	0.24	0.18	[0.01, 0.64]	NA	NA	NA	NA	NA	1.000	15429	19416
sd(Daily individual’s experienced pressure)	0.74	0.42	[0.07, 1.75]	NA	NA	NA	NA	NA	1.000	12028	12872
sd(Daily partner’s experienced pressure)	0.75	0.61	[0.04, 2.34]	NA	NA	NA	NA	NA	1.000	11524	18542
sd(Daily individual’s experienced pushing)	0.26	0.22	[0.01, 0.79]	NA	NA	NA	NA	NA	1.000	11607	19536
sd(Daily partner’s experienced pushing)	0.20	0.18	[0.01, 0.75]	NA	NA	NA	NA	NA	1.000	18835	21947
sd(barriers_self_cw)	0.75	0.52	[0.04, 1.86]	NA	NA	NA	NA	NA	1.000	10204	16342
sd(facilitators_self_cw)	1.15	0.44	[0.23, 2.08]	NA	NA	NA	NA	NA	1.001	9945	7825
sd(persuasion_self_cw:barriers_self_cw)	0.37	0.25	[0.02, 0.93]	NA	NA	NA	NA	NA	1.000	13473	16706
sd(barriers_self_cw:pressure_self_cw)	0.41	0.37	[0.02, 1.62]	NA	NA	NA	NA	NA	1.000	21785	20837
sd(barriers_self_cw:pushing_self_cw)	0.30	0.25	[0.01, 0.94]	NA	NA	NA	NA	NA	1.000	14225	20193
sd(persuasion_self_cw:facilitators_self_cw)	0.21	0.18	[0.01, 0.70]	NA	NA	NA	NA	NA	1.001	14719	19964
sd(pressure_self_cw:facilitators_self_cw)	1.19	0.62	[0.16, 2.69]	NA	NA	NA	NA	NA	1.000	13974	13085
sd(pushing_self_cw:facilitators_self_cw)	0.35	0.29	[0.02, 1.12]	NA	NA	NA	NA	NA	1.000	13347	21462
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
disc	1.00	0.00	[1.00, 1.00]	NA	NA	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()

## Possible multicollinearity between b_Intercept[4] and b_Intercept[2]
##   (r = 0.79), b_Intercept[4] and b_Intercept[3] (r = 0.85),
##   b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.76). This
##   might lead to inappropriate results. See 'Details' in '?rope'.

#conditional_spaghetti(
#  reactance_ordinal, 
#  effects = c('persuasion_self_cw', 'pressure_self_cw')
#  , group_var = 'coupleID'
#  #, n_groups = 15
#  , plot_full_range = T
#)

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cb + 
    facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw * barriers_self_cw + 
    pressure_self_cw * barriers_self_cw + 
    pushing_self_cw * barriers_self_cw + 
    
    persuasion_self_cw * facilitators_self_cw + 
    pressure_self_cw * facilitators_self_cw + 
    pushing_self_cw * facilitators_self_cw + 
    
    persuasion_partner_cw +
    pressure_partner_cw +
    pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

is_reactance_digest <- digest::digest(is_reactance)

if (check_models) {
  check_brms(is_reactance)
  DHARMa.check_brms.all(is_reactance, integer = FALSE)
}

## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.12 [1.08, 1.19]         1.06      0.89
##       barriers_self_cw 1.56 [1.48, 1.66]         1.25      0.64
##       pressure_self_cw 1.06 [1.02, 1.13]         1.03      0.95
##        pushing_self_cw 1.13 [1.09, 1.19]         1.06      0.89
##   facilitators_self_cw 1.81 [1.71, 1.93]         1.35      0.55
##  persuasion_partner_cw 1.28 [1.22, 1.35]         1.13      0.78
##    pressure_partner_cw 1.06 [1.02, 1.13]         1.03      0.95
##     pushing_partner_cw 1.24 [1.19, 1.31]         1.11      0.81
##     persuasion_self_cb 3.74 [3.49, 4.03]         1.93      0.27
##  persuasion_partner_cb 4.03 [3.75, 4.34]         2.01      0.25
##       pressure_self_cb 2.27 [2.13, 2.43]         1.51      0.44
##    pressure_partner_cb 2.22 [2.09, 2.38]         1.49      0.45
##        pushing_self_cb 2.36 [2.21, 2.52]         1.54      0.42
##     pushing_partner_cb 2.43 [2.28, 2.60]         1.56      0.41
##       barriers_self_cb 1.26 [1.20, 1.33]         1.12      0.80
##   facilitators_self_cb 1.14 [1.10, 1.21]         1.07      0.88
##                    day 1.07 [1.04, 1.14]         1.04      0.93
##  Tolerance 95% CI
##      [0.84, 0.92]
##      [0.60, 0.67]
##      [0.89, 0.98]
##      [0.84, 0.92]
##      [0.52, 0.59]
##      [0.74, 0.82]
##      [0.89, 0.98]
##      [0.76, 0.84]
##      [0.25, 0.29]
##      [0.23, 0.27]
##      [0.41, 0.47]
##      [0.42, 0.48]
##      [0.40, 0.45]
##      [0.38, 0.44]
##      [0.75, 0.83]
##      [0.83, 0.91]
##      [0.88, 0.96]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        normal         31%
##          beta         19%
##        cauchy         19%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##      bernoulli         75%
##  beta-binomial         16%
##       binomial          9%
## 
## Divergences:

## 0 of 40000 iterations ended with a divergence.

## 
## Tree depth:

## 0 of 40000 iterations saturated the maximum tree depth of 10.

## 
## Energy:

## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 115 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 486, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.000000000 0.007561553
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                                    0

if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}

summary_is_reactance <- summarize_brms(
  is_reactance, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 

## Sampling priors, please wait...

summary_is_reactance %>%
  print_df(rows_to_pack = rows_to_pack)

	OR	SE	95% CI	pd	ROPE	inside ROPE	BF	BF_Evidence	Rhat	Bulk_ESS	Tail_ESS
Intercept	0.40	0.23	[0.12, 1.25]	0.944	[0.83, 1.20]	0.072	0.353	Weak Evidence for Null	1.000	19541	24420
Within-Person Effects
Daily individual’s experienced persuasion	0.60*	0.14	[0.36, 0.93]	0.987	[0.83, 1.20]	0.066	1.694	Weak Evidence	1.000	18821	21645
Daily partner’s experienced persuasion	1.35	0.36	[0.82, 2.51]	0.877	[0.83, 1.20]	0.295	0.225	Moderate Evidence for Null	1.000	15298	17996
Daily individual’s experienced pressure	6.78*	6.14	[1.40, 63.18]	0.991	[0.83, 1.20]	0.011	3.813	Moderate Evidence	1.000	13931	18160
Daily partner’s experienced pressure	1.51	1.23	[0.29, 13.36]	0.706	[0.83, 1.20]	0.162	0.205	Moderate Evidence for Null	1.000	16765	18592
Daily individual’s experienced pushing	1.59	0.41	[0.99, 2.98]	0.973	[0.83, 1.20]	0.118	0.678	Weak Evidence for Null	1.000	17871	19024
Daily partner’s experienced pushing	0.87	0.28	[0.45, 1.73]	0.670	[0.83, 1.20]	0.391	0.137	Moderate Evidence for Null	1.000	21742	23601
Day	1.65	1.26	[0.37, 7.62]	0.745	[0.83, 1.20]	0.153	0.103	Moderate Evidence for Null	1.000	26409	28870
Own actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Partner actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily barriers	1.02	0.56	[0.34, 3.35]	0.518	[0.83, 1.20]	0.259	0.096	Strong Evidence for Null	1.000	25169	27753
Daily facilitators	1.08	0.64	[0.30, 3.40]	0.549	[0.83, 1.20]	0.237	0.115	Moderate Evidence for Null	1.000	20365	24604
Daily wear time	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean individual’s experienced persuasion	4.44	6.70	[0.24, 116.32]	0.845	[0.83, 1.20]	0.062	0.256	Moderate Evidence for Null	1.000	12271	19664
Mean partner’s experienced persuasion	8.48	13.61	[0.43, 277.62]	0.918	[0.83, 1.20]	0.037	0.371	Weak Evidence for Null	1.000	11662	18399
Mean individual’s experienced pressure	59.85	151.80	[0.54, 15556.44]	0.955	[0.83, 1.20]	0.015	0.708	Weak Evidence for Null	1.000	14113	19622
Mean partner’s experienced pressure	0.73	2.04	[0.00, 144.08]	0.545	[0.83, 1.20]	0.051	0.209	Moderate Evidence for Null	1.000	9900	18663
Mean individual’s experienced pushing	1.59	3.71	[0.02, 178.19]	0.578	[0.83, 1.20]	0.060	0.161	Moderate Evidence for Null	1.000	11104	18225
Mean partner’s experienced pushing	0.01	0.02	[0.00, 1.03]	0.974	[0.83, 1.20]	0.008	1.108	Weak Evidence	1.000	11127	16434
Mean barriers	0.67	0.73	[0.07, 5.64]	0.645	[0.83, 1.20]	0.130	0.129	Moderate Evidence for Null	1.000	18534	26202
Mean facilitators	0.76	0.54	[0.18, 3.15]	0.648	[0.83, 1.20]	0.191	0.174	Moderate Evidence for Null	1.000	14740	22473
Mean wear time	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Interactions
persuasion_self_cw:barriers_self_cw	0.76	0.24	[0.38, 1.39]	0.821	[0.83, 1.20]	0.310	0.177	Moderate Evidence for Null	1.000	23607	22535
barriers_self_cw:pressure_self_cw	1.58	1.52	[0.22, 13.05]	0.686	[0.83, 1.20]	0.136	0.175	Moderate Evidence for Null	1.000	23408	28063
barriers_self_cw:pushing_self_cw	0.65	0.25	[0.27, 1.37]	0.880	[0.83, 1.20]	0.194	0.230	Moderate Evidence for Null	1.000	19356	22441
persuasion_self_cw:facilitators_self_cw	0.88	0.25	[0.48, 1.56]	0.679	[0.83, 1.20]	0.436	0.135	Moderate Evidence for Null	1.000	17186	19441
pressure_self_cw:facilitators_self_cw	2.55	3.26	[0.22, 45.77]	0.779	[0.83, 1.20]	0.090	0.321	Weak Evidence for Null	1.000	12990	17446
pushing_self_cw:facilitators_self_cw	0.73	0.31	[0.32, 2.07]	0.758	[0.83, 1.20]	0.232	0.246	Moderate Evidence for Null	1.000	12633	14738
Random Effects
sd(Intercept)	2.64	0.61	[1.64, 4.08]	NA	NA	NA	NA	NA	1.001	7903	16817
sd(Daily individual’s experienced persuasion)	0.66	0.36	[0.06, 1.45]	NA	NA	NA	NA	NA	1.000	4439	6807
sd(Daily partner’s experienced persuasion)	0.68	0.35	[0.09, 1.54]	NA	NA	NA	NA	NA	1.000	7782	7552
sd(Daily individual’s experienced pressure)	2.38	1.02	[0.61, 4.75]	NA	NA	NA	NA	NA	1.000	8978	7015
sd(Daily partner’s experienced pressure)	1.83	1.17	[0.14, 4.52]	NA	NA	NA	NA	NA	1.000	11193	12241
sd(Daily individual’s experienced pushing)	0.55	0.40	[0.03, 1.53]	NA	NA	NA	NA	NA	1.001	6460	10411
sd(Daily partner’s experienced pushing)	0.45	0.37	[0.02, 1.45]	NA	NA	NA	NA	NA	1.000	11809	15703
sd(barriers_self_cw)	1.13	0.71	[0.08, 2.66]	NA	NA	NA	NA	NA	1.001	7367	9451
sd(facilitators_self_cw)	1.39	0.60	[0.22, 2.73]	NA	NA	NA	NA	NA	1.000	7616	7172
sd(persuasion_self_cw:barriers_self_cw)	0.58	0.42	[0.03, 1.61]	NA	NA	NA	NA	NA	1.001	8107	12898
sd(barriers_self_cw:pressure_self_cw)	0.97	0.85	[0.04, 3.40]	NA	NA	NA	NA	NA	1.000	14182	15002
sd(barriers_self_cw:pushing_self_cw)	0.84	0.46	[0.08, 1.94]	NA	NA	NA	NA	NA	1.001	8459	9366
sd(persuasion_self_cw:facilitators_self_cw)	0.33	0.29	[0.02, 1.15]	NA	NA	NA	NA	NA	1.000	9640	17208
sd(pressure_self_cw:facilitators_self_cw)	2.54	1.21	[0.34, 5.20]	NA	NA	NA	NA	NA	1.000	10620	8039
sd(pushing_self_cw:facilitators_self_cw)	0.87	0.60	[0.06, 2.36]	NA	NA	NA	NA	NA	1.001	6284	10446
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.73). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="05_SensitivityBarriersInteraction-ONLY_DAYS_WITH_PLAN-_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw

hypothesis(is_reactance, "exp(pressure_self_cw) > exp(pushing_self_cw)")

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (exp(pressure_sel... > 0    11.61        32    -0.03    39.96      18.29
##   Post.Prob Star
## 1      0.95     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

summary_all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  is_reactance,
  
  stats_to_report = c('CI'),
  model_rows_random = model_rows_random,
  model_rows_fixed = model_rows_fixed,
  model_rownames_random = model_rownames_random,
  model_rownames_fixed = model_rownames_fixed
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate,
## stats_to_report = stats_to_report, : Coefficients were exponentiated.
## Double check if this was intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate,
## stats_to_report = stats_to_report, : Coefficients were exponentiated.
## Double check if this was intended.

[1] “mood_gauss” [1] “is_reactance”

summary_all_models <- summary_all_models %>%
  print_df(rows_to_pack = rows_to_pack) %>%
  add_header_above(
    c(
      " ", "Hurdle Component" = 2, "Non-Zero Component" = 2,
      " " = 6
    )
  ) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 4,  
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      #"Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2
    )
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack,
  file.path("Output", paste0("AllModels", suffix, ".xlsx")), 
  merge_option = 'both', 
  simplify_2nd_row = FALSE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
summary_all_models

	Subjective MVPA Hurdle Lognormal				Device-Based MVPA Log (Gaussian)		Mood Gaussian		Reactance Dichotome
	Hurdle Component		Non-Zero Component
	exp(Est.)_hu pa_sub	95% CI_hu pa_sub	exp(Est.)_nonzero pa_sub	95% CI_nonzero pa_sub	exp(Est.) pa_obj_log	95% CI pa_obj_log	Est. mood_gauss	95% CI mood_gauss	OR is_reactance	95% CI is_reactance
Intercept	3.41***	[ 1.79, 6.47]	50.01***	[42.13, 59.09]	123.06***	[108.16, 139.27]	3.79***	[ 3.58, 4.01]	0.40	[0.12, 1.25]
Within-Person Effects
Daily individual’s experienced persuasion	1.87***	[ 1.41, 2.61]	1.00	[ 0.95, 1.06]	1.02	[ 0.98, 1.06]	0.01	[-0.04, 0.06]	0.60*	[0.36, 0.93]
Daily partner’s experienced persuasion	1.40**	[ 1.12, 1.83]	1.02	[ 0.98, 1.07]	1.03	[ 0.99, 1.07]	0.01	[-0.04, 0.06]	1.35	[0.82, 2.51]
Daily individual’s experienced pressure	0.86	[ 0.39, 2.09]	0.91	[ 0.78, 1.05]	0.97	[ 0.87, 1.07]	-0.02	[-0.16, 0.12]	6.78*	[1.40, 63.18]
Daily partner’s experienced pressure	1.69	[ 0.85, 5.72]	0.97	[ 0.88, 1.06]	0.97	[ 0.90, 1.05]	-0.05	[-0.18, 0.07]	1.51	[0.29, 13.36]
Daily individual’s experienced pushing	1.24	[ 0.77, 2.01]	0.99	[ 0.93, 1.06]	1.02	[ 0.96, 1.08]	0.02	[-0.06, 0.10]	1.59	[0.99, 2.98]
Daily partner’s experienced pushing	1.60**	[ 1.13, 2.54]	0.96	[ 0.90, 1.02]	1.00	[ 0.96, 1.05]	0.09*	[ 0.02, 0.16]	0.87	[0.45, 1.73]
Day	0.91	[ 0.55, 1.51]	0.95	[ 0.83, 1.09]	0.96	[ 0.87, 1.06]	0.26***	[ 0.11, 0.41]	1.65	[0.37, 7.62]
Own actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Partner actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily barriers	1.14	[ 0.64, 2.04]	1.26***	[ 1.13, 1.41]	1.01	[ 0.94, 1.07]	-0.12	[-0.25, 0.01]	1.02	[0.34, 3.35]
Daily facilitators	1.17	[ 0.68, 1.96]	1.19**	[ 1.06, 1.33]	1.04	[ 0.99, 1.10]	0.11**	[ 0.05, 0.17]	1.08	[0.30, 3.40]
Daily wear time	NA	NA	NA	NA	1.00***	[ 1.00, 1.00]	NA	NA	NA	NA
Between-Person Effects
Mean individual’s experienced persuasion	1.07	[ 0.27, 4.02]	1.06	[ 0.75, 1.48]	0.99	[ 0.69, 1.40]	0.13	[-0.47, 0.72]	4.44	[0.24, 116.32]
Mean partner’s experienced persuasion	1.43	[ 0.38, 5.24]	0.96	[ 0.69, 1.32]	0.86	[ 0.61, 1.23]	0.42	[-0.18, 1.00]	8.48	[0.43, 277.62]
Mean individual’s experienced pressure	0.18	[ 0.03, 1.11]	0.68	[ 0.29, 1.58]	1.12	[ 0.77, 1.64]	-0.07	[-0.67, 0.53]	59.85	[0.54, 15556.44]
Mean partner’s experienced pressure	0.19	[ 0.03, 1.17]	0.49	[ 0.21, 1.15]	1.04	[ 0.72, 1.49]	-0.32	[-0.91, 0.28]	0.73	[0.00, 144.08]
Mean individual’s experienced pushing	1.32	[ 0.20, 8.70]	1.53	[ 0.84, 2.80]	1.11	[ 0.69, 1.82]	0.23	[-0.61, 1.08]	1.59	[0.02, 178.19]
Mean partner’s experienced pushing	1.45	[ 0.22, 9.51]	1.65	[ 0.91, 3.01]	1.27	[ 0.78, 2.05]	-0.06	[-0.90, 0.79]	0.01	[0.00, 1.03]
Mean barriers	0.50**	[ 0.30, 0.83]	0.97	[ 0.83, 1.12]	0.80***	[ 0.72, 0.89]	-0.50***	[-0.65, -0.34]	0.67	[0.07, 5.64]
Mean facilitators	1.24	[ 0.86, 1.80]	1.08	[ 0.98, 1.19]	0.96	[ 0.89, 1.03]	0.29***	[ 0.18, 0.39]	0.76	[0.18, 3.15]
Mean wear time	NA	NA	NA	NA	1.00*	[ 1.00, 1.00]	NA	NA	NA	NA
Interactions
persuasion_self_cw:barriers_self_cw	1.02	[ 0.65, 1.63]	1.03	[ 0.96, 1.11]	0.97	[ 0.93, 1.02]	0.00	[-0.08, 0.07]	0.76	[0.38, 1.39]
barriers_self_cw:pressure_self_cw	0.34	[ 0.09, 1.21]	0.95	[ 0.75, 1.20]	0.96	[ 0.81, 1.14]	-0.14	[-0.47, 0.20]	1.58	[0.22, 13.05]
barriers_self_cw:pushing_self_cw	2.50*	[ 1.14, 6.03]	0.94	[ 0.85, 1.05]	1.02	[ 0.95, 1.09]	0.08	[-0.02, 0.18]	0.65	[0.27, 1.37]
persuasion_self_cw:facilitators_self_cw	1.29	[ 0.87, 2.10]	1.02	[ 0.96, 1.09]	0.99	[ 0.95, 1.04]	0.00	[-0.06, 0.06]	0.88	[0.48, 1.56]
pressure_self_cw:facilitators_self_cw	0.88	[ 0.18, 4.87]	0.94	[ 0.73, 1.18]	0.92	[ 0.69, 1.14]	0.00	[-0.24, 0.31]	2.55	[0.22, 45.77]
pushing_self_cw:facilitators_self_cw	1.54	[ 0.77, 3.92]	0.95	[ 0.88, 1.04]	1.03	[ 0.96, 1.09]	0.02	[-0.08, 0.13]	0.73	[0.32, 2.07]
Random Effects
sd(Intercept)	1.64	[1.23, 2.23]	0.32	[0.25, 0.43]	0.34	[0.26, 0.45]	0.62	[0.49, 0.81]	2.64	[1.64, 4.08]
sd(Daily individual’s experienced persuasion)	0.21	[0.01, 0.61]	0.08	[0.03, 0.14]	0.05	[0.01, 0.09]	0.04	[0.00, 0.11]	0.66	[0.06, 1.45]
sd(Daily partner’s experienced persuasion)	0.30	[0.03, 0.66]	0.06	[0.01, 0.11]	0.05	[0.01, 0.10]	0.04	[0.00, 0.11]	0.68	[0.09, 1.54]
sd(Daily individual’s experienced pressure)	0.54	[0.02, 2.15]	0.08	[0.00, 0.33]	0.05	[0.00, 0.20]	0.07	[0.00, 0.25]	2.38	[0.61, 4.75]
sd(Daily partner’s experienced pressure)	0.77	[0.04, 2.37]	0.05	[0.00, 0.18]	0.03	[0.00, 0.12]	0.09	[0.00, 0.29]	1.83	[0.14, 4.52]
sd(Daily individual’s experienced pushing)	0.70	[0.07, 1.53]	0.05	[0.00, 0.14]	0.08	[0.01, 0.15]	0.09	[0.01, 0.20]	0.55	[0.03, 1.53]
sd(Daily partner’s experienced pushing)	0.44	[0.03, 1.23]	0.06	[0.00, 0.14]	0.04	[0.00, 0.11]	0.09	[0.01, 0.19]	0.45	[0.02, 1.45]
sd(barriers_self_cw)	1.51	[1.10, 2.08]	0.23	[0.14, 0.35]	0.13	[0.06, 0.20]	0.32	[0.23, 0.43]	1.13	[0.08, 2.66]
sd(facilitators_self_cw)	1.30	[0.84, 1.95]	0.27	[0.19, 0.39]	0.07	[0.01, 0.15]	0.05	[0.00, 0.16]	1.39	[0.22, 2.73]
sd(persuasion_self_cw:barriers_self_cw)	0.58	[0.05, 1.26]	0.07	[0.00, 0.18]	0.03	[0.00, 0.09]	0.07	[0.00, 0.18]	0.58	[0.03, 1.61]
sd(barriers_self_cw:pressure_self_cw)	0.67	[0.03, 2.58]	0.11	[0.01, 0.42]	0.07	[0.00, 0.31]	0.41	[0.08, 0.85]	0.97	[0.04, 3.40]
sd(barriers_self_cw:pushing_self_cw)	1.26	[0.38, 2.39]	0.15	[0.05, 0.26]	0.03	[0.00, 0.10]	0.05	[0.00, 0.18]	0.84	[0.08, 1.94]
sd(persuasion_self_cw:facilitators_self_cw)	0.59	[0.11, 1.23]	0.07	[0.00, 0.16]	0.03	[0.00, 0.09]	0.04	[0.00, 0.13]	0.33	[0.02, 1.15]
sd(pressure_self_cw:facilitators_self_cw)	1.60	[0.10, 4.28]	0.12	[0.01, 0.51]	0.24	[0.02, 0.62]	0.16	[0.01, 0.55]	2.54	[0.34, 5.20]
sd(pushing_self_cw:facilitators_self_cw)	0.94	[0.08, 2.45]	0.05	[0.00, 0.17]	0.04	[0.00, 0.14]	0.11	[0.01, 0.26]	0.87	[0.06, 2.36]
Additional Parameters
sigma	NA	NA	0.62	[0.59, 0.65]	0.53	[0.51, 0.55]	0.84	[0.82, 0.88]	NA	NA

Plot Interaction effect.

1. Probability of being active (crossing the hurdle)

library(ggplot2)
library(plotly)


p <- conditional_effects(
  pa_sub, 
  dpar = 'hu',
  effects = 'pushing_self_cw:barriers_self_cw',
  re_formula = NA, 
  method = 'posterior_linpred',
  plot = FALSE
)

p2 <- p

p2[[1]]$estimate__ <- 1 - p[[1]]$estimate__
p2[[1]]$lower__ <- 1 - p[[1]]$upper__
p2[[1]]$upper__ <- 1 - p[[1]]$lower__

plots <- plot(p2, theme = theme_minimal(), plot = FALSE) 

plots[[1]] + 
  scale_color_manual(values = c("steelblue", "forestgreen", "coral")) +
  scale_fill_manual(values = c("steelblue", "forestgreen", "coral"))

Test slope of pushing at high levels of barriers

library(tidybayes)
library(dplyr)
library(ggplot2)

plot_slope_by_moderator <- function(
  mod,            
  df,             
  moderator,      
  predictor_par,  
  interaction_par,
  exponentiate = FALSE,   
  n_points = 1000,
  fill_colors = c("Significant" = "steelblue",
                  "Non-significant" = "coral"),
  line_colors = fill_colors,  # Reuse fill colors by default
  vjust_val = -1.5
) {
  
  # Generate a sequence of moderator values
  moderator_seq <- seq(
    from = min(df[[moderator]], na.rm = TRUE),
    to = max(df[[moderator]], na.rm = TRUE),
    length.out = n_points
  )
  
  # Extract posterior draws for the parameters of interest
  library(rlang)
  predictor_sym <- rlang::sym(predictor_par)
  interaction_sym <- rlang::sym(interaction_par)

  post_draws <- mod %>%
    spread_draws(!!predictor_sym, !!interaction_sym)
  
  # For each moderator value, compute the slope for each posterior draw
  slope_list <- lapply(moderator_seq, function(mod_val) {
    tibble(
      moderator_value = mod_val,
      slope = - (post_draws[[predictor_par]] + mod_val * post_draws[[interaction_par]])
    )
  })
  
  # Combine the list into a single data frame
  slope_df <- bind_rows(slope_list)
  
  # Summarize the posterior for each moderator value (using a 95% credible interval)
  slope_summary <- slope_df %>%
    group_by(moderator_value) %>%
    summarize(
      slope_median = median(slope),
      slope_lower  = quantile(slope, probs = 0.025),
      slope_upper  = quantile(slope, probs = 0.975)
    ) %>%
    # Create significance flag: "sig" if CI excludes 0, "ns" otherwise
    mutate(sig_flag = if_else(slope_lower > 0 | slope_upper < 0, "sig", "ns")) %>%
    arrange(moderator_value)
  
  # If exponentiate == TRUE, transform slope columns
  if (exponentiate) {
    slope_summary <- slope_summary %>%
      mutate(
        slope_median_exp = exp(slope_median),
        slope_lower_exp  = exp(slope_lower),
        slope_upper_exp  = exp(slope_upper)
      )
  }
  
  # Identify boundary rows and duplicate them for each category
  slope_summary <- slope_summary %>%
    mutate(
      next_sig_flag = lead(sig_flag),
      is_boundary   = sig_flag != next_sig_flag & !is.na(next_sig_flag)
    )
  
  boundary_rows <- slope_summary %>%
    filter(is_boundary)
  
  boundary_rows_next <- boundary_rows
  boundary_rows_next$sig_flag <- boundary_rows_next$next_sig_flag
  
  slope_summary_extended <- bind_rows(slope_summary, boundary_rows_next)
  
  # Split into "sig" and "ns"
  slope_summary_sig <- slope_summary_extended %>% filter(sig_flag == "sig")
  slope_summary_ns  <- slope_summary_extended %>% filter(sig_flag == "ns")
  
  # Build the base plot:
  #    We'll use either slope_{median, lower, upper} or slope_{median_exp, lower_exp, upper_exp} 
  #    depending on exponentiate.
  if (exponentiate) {
    median_col <- "slope_median_exp"
    lower_col <- "slope_lower_exp"
    upper_col <- "slope_upper_exp"
    y_lab <- "Odds Ratio (Exponentiated Slope)"
    ref_line <- 1  # reference line for ORs is 1
  } else {
    median_col <- "slope_median"
    lower_col <- "slope_lower"
    upper_col <- "slope_upper"
    y_lab <- "Slope (Log-Odds Scale)"
    ref_line <- 0
  }
  
  p <- ggplot() +
    # Non-significant region
    geom_ribbon(
      data = slope_summary_ns,
      aes_string(
        x = "moderator_value",
        ymin = lower_col,
        ymax = upper_col,
        fill = '"Non-significant"'
      ),
      alpha = 0.2
    ) +
    geom_line(
      data = slope_summary_ns,
      aes_string(
        x = "moderator_value",
        y = median_col,
        color = '"Non-significant"'
      ),
      size = 1
    ) +
    # Significant region
    geom_ribbon(
      data = slope_summary_sig,
      aes_string(
        x = "moderator_value",
        ymin = lower_col,
        ymax = upper_col,
        fill = '"Significant"'
      ),
      alpha = 0.2
    ) +
    geom_line(
      data = slope_summary_sig,
      aes_string(
        x = "moderator_value",
        y = median_col,
        color = '"Significant"'
      ),
      size = 1
    ) +
    # Reference line (0 for log-odds, 1 for odds-ratios)
    geom_hline(
      yintercept = ref_line, 
      linetype = "dashed"
    ) +
    # Manual scales for fill & color
    scale_fill_manual(
      values = fill_colors
    ) +
    scale_color_manual(
      values = line_colors
    ) +
    labs(
      x = moderator,
      y = y_lab,
      fill = "Significance",
      color = "Significance"
    ) +
    theme_minimal()
  
  # Add boundary lines (where significance changes) & labels
  p2 <- p +
    geom_vline(
      data = boundary_rows,
      aes(xintercept = moderator_value),
      color = fill_colors["Significant"],
      alpha = 0.35,
      size = 0.7,
      linetype = "dashed"
    ) +
    geom_text(
      data = boundary_rows,
      aes(x = moderator_value, y = -Inf, label = round(moderator_value, 3)),
      vjust = vjust_val,
      size  = 3.7
    )
  
  return(p2)
}


p_log_odds <- plot_slope_by_moderator(
  mod = pa_sub, 
  df = df_double, 
  moderator = "barriers_self_cw", 
  predictor_par = "b_hu_pushing_self_cw", 
  interaction_par = "b_hu_barriers_self_cw:pushing_self_cw",
  exponentiate = FALSE
)

p_log_odds

p_odds_ratio <- plot_slope_by_moderator(
  mod = pa_sub, 
  df = df_double, 
  moderator = "barriers_self_cw", 
  predictor_par = "b_hu_pushing_self_cw", 
  interaction_par = "b_hu_barriers_self_cw:pushing_self_cw",
  exponentiate = TRUE,
  vjust_val = -12,
) + coord_cartesian(ylim = c(1, 18))

p_odds_ratio

# Export
ggsave(
  filename = file.path("Output", "InteractPlot_log_odds.png"), 
  plot = p_log_odds,
  width = 8,
  height = 5,
  dpi = 300,
  bg = 'white'
)

ggsave(
  filename = file.path("Output", "InteractPlot_odds_ratio.png"), 
  plot = p_odds_ratio, 
  width = 8, 
  height = 5, 
  dpi = 300,
  bg = 'white'
)

2. Expected MVPA when the hurdle was crossed (interaction NOT significant)

p3 <- conditional_effects(
  pa_sub, 
  dpar = 'mu',
  effects = 'pushing_self_cw:barriers_self_cw',
  #re_formula = NULL, 
  plot = FALSE
)

p3[[1]]$estimate__ <- exp(p3[[1]]$estimate__)
p3[[1]]$lower__ <- exp(p3[[1]]$lower__)
p3[[1]]$upper__ <- exp(p3[[1]]$upper__)

plot(p3, ask = FALSE)

2. Expected physical activity in Minutes (when both components combined)

# Expected physical activity
p4 <- conditional_effects(
  pa_sub, 
  effects = 'pushing_self_cw:barriers_self_cw',
  #re_formula = NULL,
  plot = FALSE
)

plot(p4, ask = FALSE)

report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)

report::report_packages()

• rmarkdown (version 2.29; Allaire J et al., 2024)
• beepr (version 2.0; Bååth R, 2024)
• R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
• R.oo (version 1.27.0; Bengtsson H, 2003)
• R.utils (version 2.12.3; Bengtsson H, 2023)
• brms (version 2.22.0; Bürkner P, 2017)
• posterior (version 1.6.0; Bürkner P et al., 2024)
• extrafont (version 0.19; Chang W, 2023)
• digest (version 0.6.37; Eddelbuettel D, 2024)
• Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
• bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
• lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
• DHARMa (version 0.4.7; Hartig F, 2024)
• rlang (version 1.1.4; Henry L, Wickham H, 2024)
• tidybayes (version 3.0.7; Kay M, 2024)
• wbCorr (version 0.1.22; Küng P, 2023)
• see (version 0.9.0; Lüdecke D et al., 2021)
• tibble (version 3.2.1; Müller K, Wickham H, 2023)
• patchwork (version 1.3.0; Pedersen T, 2024)
• R (version 4.4.2; R Core Team, 2024)
• openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
• plotly (version 4.10.4; Sievert C, 2020)
• ggplot2 (version 3.5.1; Wickham H, 2016)
• forcats (version 1.0.0; Wickham H, 2023)
• stringr (version 1.5.1; Wickham H, 2023)
• rvest (version 1.0.4; Wickham H, 2024)
• tidyverse (version 2.0.0; Wickham H et al., 2019)
• readxl (version 1.4.3; Wickham H, Bryan J, 2023)
• dplyr (version 1.1.4; Wickham H et al., 2023)
• purrr (version 1.0.2; Wickham H, Henry L, 2023)
• readr (version 2.1.5; Wickham H et al., 2024)
• xml2 (version 1.3.6; Wickham H et al., 2023)
• tidyr (version 1.3.1; Wickham H et al., 2024)
• ggridges (version 0.5.6; Wilke C, 2024)
• knitr (version 1.49; Xie Y, 2024)
• kableExtra (version 1.4.0; Zhu H, 2024)

report::cite_packages()

• Allaire J, Xie Y, Dervieux C, McPherson J, Luraschi J, Ushey K, Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2024). rmarkdown: Dynamic Documents for R. R package version 2.29, https://github.com/rstudio/rmarkdown. Xie Y, Allaire J, Grolemund G (2018). R Markdown: The Definitive Guide. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9781138359338, https://bookdown.org/yihui/rmarkdown. Xie Y, Dervieux C, Riederer E (2020). R Markdown Cookbook. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9780367563837, https://bookdown.org/yihui/rmarkdown-cookbook.
• Bååth R (2024). beepr: Easily Play Notification Sounds on any Platform. R package version 2.0, https://github.com/rasmusab/beepr.
• Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
• Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
• Bengtsson H (2023). R.utils: Various Programming Utilities. R package version 2.12.3, https://github.com/HenrikBengtsson/R.utils, https://henrikbengtsson.github.io/R.utils/.
• Bürkner P (2017). “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software, 80(1), 1-28. doi:10.18637/jss.v080.i01 https://doi.org/10.18637/jss.v080.i01. Bürkner P (2018). “Advanced Bayesian Multilevel Modeling with the R Package brms.” The R Journal, 10(1), 395-411. doi:10.32614/RJ-2018-017 https://doi.org/10.32614/RJ-2018-017. Bürkner P (2021). “Bayesian Item Response Modeling in R with brms and Stan.” Journal of Statistical Software, 100(5), 1-54. doi:10.18637/jss.v100.i05 https://doi.org/10.18637/jss.v100.i05.
• Bürkner P, Gabry J, Kay M, Vehtari A (2024). “posterior: Tools for Working with Posterior Distributions.” R package version 1.6.0, https://mc-stan.org/posterior/. Vehtari A, Gelman A, Simpson D, Carpenter B, Bürkner P (2021). “Rank-normalization, folding, and localization: An improved Rhat for assessing convergence of MCMC (with discussion).” Bayesian Analysis.
• Chang W (2023). extrafont: Tools for Using Fonts. R package version 0.19, https://github.com/wch/extrafont.
• Eddelbuettel D (2024). digest: Create Compact Hash Digests of R Objects. R package version 0.6.37, https://dirk.eddelbuettel.com/code/digest.html, https://github.com/eddelbuettel/digest.
• Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2024). Rcpp: Seamless R and C++ Integration. R package version 1.0.13-1, https://dirk.eddelbuettel.com/code/rcpp.html, https://github.com/RcppCore/Rcpp, https://www.rcpp.org. Eddelbuettel D, François R (2011). “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software, 40(8), 1-18. doi:10.18637/jss.v040.i08 https://doi.org/10.18637/jss.v040.i08. Eddelbuettel D (2013). Seamless R and C++ Integration with Rcpp. Springer, New York. doi:10.1007/978-1-4614-6868-4 https://doi.org/10.1007/978-1-4614-6868-4, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). “Extending R with C++: A Brief Introduction to Rcpp.” The American Statistician, 72(1), 28-36. doi:10.1080/00031305.2017.1375990 https://doi.org/10.1080/00031305.2017.1375990.
• Gabry J, Mahr T (2024). “bayesplot: Plotting for Bayesian Models.” R package version 1.11.1, https://mc-stan.org/bayesplot/. Gabry J, Simpson D, Vehtari A, Betancourt M, Gelman A (2019). “Visualization in Bayesian workflow.” J. R. Stat. Soc. A, 182, 389-402. doi:10.1111/rssa.12378 https://doi.org/10.1111/rssa.12378.
• Grolemund G, Wickham H (2011). “Dates and Times Made Easy with lubridate.” Journal of Statistical Software, 40(3), 1-25. https://www.jstatsoft.org/v40/i03/.
• Hartig F (2024). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models. R package version 0.4.7, http://florianhartig.github.io/DHARMa/.
• Henry L, Wickham H (2024). rlang: Functions for Base Types and Core R and ‘Tidyverse’ Features. R package version 1.1.4, https://github.com/r-lib/rlang, https://rlang.r-lib.org.
• Kay M (2024). tidybayes: Tidy Data and Geoms for Bayesian Models. doi:10.5281/zenodo.1308151 https://doi.org/10.5281/zenodo.1308151, R package version 3.0.7, http://mjskay.github.io/tidybayes/.
• Küng P (2023). wbCorr: Bivariate Within- and Between-Cluster Correlations. University of Zürich. R package version 0.1.22. https://github.com/Pascal-Kueng/wbCorr.
• Lüdecke D, Patil I, Ben-Shachar M, Wiernik B, Waggoner P, Makowski D (2021). “see: An R Package for Visualizing Statistical Models.” Journal of Open Source Software, 6(64), 3393. doi:10.21105/joss.03393 https://doi.org/10.21105/joss.03393.
• Müller K, Wickham H (2023). tibble: Simple Data Frames. R package version 3.2.1, https://github.com/tidyverse/tibble, https://tibble.tidyverse.org/.
• Pedersen T (2024). patchwork: The Composer of Plots. R package version 1.3.0, https://github.com/thomasp85/patchwork, https://patchwork.data-imaginist.com.
• R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
• Schauberger P, Walker A (2024). openxlsx: Read, Write and Edit xlsx Files. R package version 4.2.7.1, https://github.com/ycphs/openxlsx, https://ycphs.github.io/openxlsx/index.html.
• Sievert C (2020). Interactive Web-Based Data Visualization with R, plotly, and shiny. Chapman and Hall/CRC. ISBN 9781138331457, https://plotly-r.com.
• Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.
• Wickham H (2023). forcats: Tools for Working with Categorical Variables (Factors). R package version 1.0.0, https://github.com/tidyverse/forcats, https://forcats.tidyverse.org/.
• Wickham H (2023). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.1, https://github.com/tidyverse/stringr, https://stringr.tidyverse.org.
• Wickham H (2024). rvest: Easily Harvest (Scrape) Web Pages. R package version 1.0.4, https://github.com/tidyverse/rvest, https://rvest.tidyverse.org/.
• Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.
• Wickham H, Bryan J (2023). readxl: Read Excel Files. R package version 1.4.3, https://github.com/tidyverse/readxl, https://readxl.tidyverse.org.
• Wickham H, François R, Henry L, Müller K, Vaughan D (2023). dplyr: A Grammar of Data Manipulation. R package version 1.1.4, https://github.com/tidyverse/dplyr, https://dplyr.tidyverse.org.
• Wickham H, Henry L (2023). purrr: Functional Programming Tools. R package version 1.0.2, https://github.com/tidyverse/purrr, https://purrr.tidyverse.org/.
• Wickham H, Hester J, Bryan J (2024). readr: Read Rectangular Text Data. R package version 2.1.5, https://github.com/tidyverse/readr, https://readr.tidyverse.org.
• Wickham H, Hester J, Ooms J (2023). xml2: Parse XML. R package version 1.3.6, https://github.com/r-lib/xml2, https://xml2.r-lib.org/.
• Wickham H, Vaughan D, Girlich M (2024). tidyr: Tidy Messy Data. R package version 1.3.1, https://github.com/tidyverse/tidyr, https://tidyr.tidyverse.org.
• Wilke C (2024). ggridges: Ridgeline Plots in ‘ggplot2’. R package version 0.5.6, https://wilkelab.org/ggridges/.
• Xie Y (2024). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.49, https://yihui.org/knitr/. Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/. Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.
• Zhu H (2024). kableExtra: Construct Complex Table with ‘kable’ and Pipe Syntax. R package version 1.4.0, https://github.com/haozhu233/kableExtra, http://haozhu233.github.io/kableExtra/.